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Abstract — This  paper  considers  the  problem  of  shooter  local¬ 
ization  using  a  network  of  soldier-worn  gunfire  detection  systems. 
Proposed  scheme  utilizes  the  benefits  of  sensor  network  layout 
of  all  the  sensors  within  a  small  combat  unit  to  help  refine 
localization  accuracy.  If  the  soldier  is  within  the  field  of  view 
of  the  shockwave,  then  using  the  acoustic  phenomena  analysis 
of  small-arms  fire,  the  gunfire  detection  system  can  localize 
the  source  of  the  incoming  fire  and  the  bullet’s  trajectory  with 
respect  to  the  sensor  location.  These  individual  solutions,  usually 
in  the  form  of  a  bearing  and  range  relative  to  the  soldier,  are 
then  relayed  to  the  central  node.  At  the  central  node  level,  the 
individual  solutions  are  fused  along  with  the  GPS  locations  on 
the  soldiers  to  yield  a  highly  accurate  geo-rectified  solution. 
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I.  Introduction 

There  is  an  eminent  need  for  highly  accurate  small-arms 
gunfire  detection  systems  on  individual  soldiers  for  added 
battlefield  situational  awareness  and  threat  assessment.  Today, 
several  acoustic  shooter  localization  systems  are  commer¬ 
cially  available  [  1 J.  [2];  an  overview  of  such  systems  can  be 
found  in  [3].  Currently  operational  Soldier  Wearable  Gunfire 
Detection  Systems  (SW-GDSs)  can  provide  an  appropriate 
level  of  localization  accuracy  as  long  as  the  soldier  is  at 
an  ideal  location  (range,  attitude,  etc.)  when  incoming  fire 
is  received  [4]-[6].  The  localization  system  suffers  severe 
performance  degradation  when  the  soldier  is  at  a  non-ideal 
location.  Moreover,  when  a  relative  solution,  i.e.,  the  shooter 
location  relative  to  the  sensor,  is  transformed  into  a  geo- 
rectified  solution  using  a  magnetometer  and  GPS,  the  solu¬ 
tion  often  becomes  unusable  due  to  localization  errors.  Geo- 
rectified  solutions  are  necessary  when  displaying  hostile  fire 
icons  on  a  Command  and  Control  Geographic  Information 
System  (C2  GIS)  map  display. 

SW-GDSs  use  acoustic  phenomena  analysis  of  small-arms 
fire  to  localize  the  source  of  incoming  fire,  usually  with  a 
bearing  and  range  relative  to  the  user  [7],  These  individual 
SW-GDSs  operate  separately  and  are  not  designed  to  exploit 
the  sensor  network  layout  of  all  the  soldiers  within  a  Small 
Combat  Unit  (SCU)  to  help  increase  accuracy.  Researchers  are 
exploring  some  novel  solutions  that  utilize  the  team  aspect  of 
these  SCUs  by  exploiting  all  SW-GDSs  in  a  squad/platoon 
to  increase  detection  rates  and  accuracy  [8] — [10].  This  paper 
presents  the  development  of  a  sensor  fusion  module  that  would 
take  full  advantage  of  the  team  aspect  of  a  SCU  to  provide  a 


fused  solution  that  would  be  highly  accurate  and  suitable  for 
a  C2  GIS  map  display  compared  to  the  individual  soldier’s 
solution.  The  objective  here  is  to  improve  accuracy  across 
an  entire  SCU  so  even  soldiers  in  non-ideal  settings  (out  of 
range,  bad  angle,  etc.)  can  exploit  the  good  solutions  from 
their  neighbors  to  come  up  with  improved  solutions:  both  geo- 
rectified  and  relative. 

The  individual  SW-GDSs  considered  here  is  composed  of 
a  passive  array  of  microphones  that  is  able  to  localize  a 
gunfire  event  by  measuring  the  direction  of  arrival  for  both 
the  acoustic  wave  generated  by  the  muzzle  blast  and  the 
shockwave  generated  by  the  supersonic  bullet  [1],  [2].  After 
detecting  a  gunfire,  the  individual  sensors  report  their  solution 
along  with  their  GPS  positions  to  a  central  node.  At  the 
central  node,  the  individual  solutions  are  fused  along  with 
the  GPS  positions  to  yield  an  highly  accurate,  geo-rectified 
solution,  which  is  then  relayed  back  to  individual  soldiers 
for  added  situational  awareness.  Structure  of  this  paper  is  as 
follows:  section  II  presents  the  measurement  model  for  the 
acoustic  sensor  nodes,  and  section  III  presents  the  localization 
algorithm  that  converts  the  sensor  measurements  to  a  gunfire 
position  estimate.  Details  of  the  central  node  data  fusion  and 
the  corresponding  nonlinear  least-squares  problem  is  given 
in  section  IV.  Section  V  presents  the  Gauss-Newton  method 
to  solve  the  nonlinear  least-squares  problem  and  section  VI 
presents  the  results  from  numerical  simulations.  Finally,  sec¬ 
tion  VII  concludes  the  paper  and  discusses  the  current  research 
challenges. 


II.  Problem  Setup 


Consider  a  SCU  consist  of  n  individual  soldiers  equipped 
with  the  SW-GDS.  In  order  to  set  up  the  problem  and  develop 
a  sensor  model,  we  first  consider  a  scenario  where  there  is 
only  one  shooter  and  the  SW-GDS  receives  both  the  muzzle 
blast  and  the  shockwave.  The  shooter  or  the  target  location  and 
the  soldier  or  the  ith  sensor  location  are  defined  as  T  and  S», 
respectively.  For  simplicity,  the  problem  is  formulated  in  SR2, 


T  e  SR2  = 


and  Si  €  SR2  = 


individual  range,  rt,  and  bearing, 
node  and  the  target  as 


.  Now  define  the 


c i>i ,  between  the  ith  sensor 


ri  =  \/(Tx-Sioc)2  +  (Ty-Siy)2  (1) 

<f>i  =  arctan  2  (Ty  -  Siy ,  Tx  -  Six)  (2) 
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Figure  1.  Geometry  of  the  bullet  trajectory  and  propagation  of  the  muzzle 
blast  and  shockwave  to  the  sensor  node. 


Figure  2.  Muzzle  blast  and  shockwave  field  of  view 


When  a  gun  fires,  the  blast  from  the  muzzle  produces  a 
spherical  acoustic  wave  that  can  be  heard  in  any  direction.  The 
bullet  travels  at  supersonic  speeds  and  produces  an  acoustic 
shockwave  that  emanates  as  a  cone  from  the  trajectory  of  the 
bullet.  Because  the  bullet  is  traveling  faster  than  the  speed 
of  sound,  the  shockwave  arrives  at  the  sensor  node  before 
the  wave  from  the  muzzle  blast,  which  we  simply  refer  to 
as  the  muzzle  blast.  Figure  1  illustrates  the  geometry  of  the 
shockwave  and  the  muzzle  blast  for  the  ith  sensor  node  when 
the  orientation  of  the  bullet  trajectory  is  u>  with  respect  to  the 
horizontal  axis.  As  the  bullet  pushes  air,  it  creates  an  impulse 
wave.  The  wavefront  is  a  cone  whose  angle  9  with  respect  to 
the  trajectory  is 


0  =  arcsin 


(3) 


where  m  is  the  Mach  number.  The  Mach  number  is  assumed 
to  be  known  since  the  typical  value  for  a  Mach  number  is 
m  =  2  [71.  Since  the  Mach  number  directly  influences  the 
range  estimates,  uncertainty  in  bullet  speed  may  be  treated  as 
range  estimation  error.  As  indicated  in  figure  1,  the  angle  fa 
indicates  the  direction  of  arrival  (DOA)  of  the  muzzle  blast, 
and  ifi  indicates  the  DOA  of  the  shockwave.  The  muzzle  blast 
DOA  is  measured  counter-clockwise  such  that  0  <  fa  <  2nl . 
For  a  more  detailed  description  of  the  scenario,  please  refer 
to  [7].  Figure  2  indicates  the  field  of  view  (FOV)  for  both  the 
muzzle  blast  and  the  shockwave.  Note  that  the  FOV  of  the 
muzzle  blast  is  27T,  i.e.,  omnidirectional,  and  the  FOV  for  the 
shockwave  is  n  —  26.  SW-GDS  receives  the  shockwave  only 
if  the  muzzle  blast  DOA  is  within  the  bounds 


7t/2  +  9  +  u>  <C  ^>2  <  37t /2  —  6  uj  (4) 

Now  the  DOA  angle  for  the  shockwave  can  be  written  as 

^  _f  -f- 6 +  0J,  if  tt  +  u>  <  fa  <  ^  —  6  +  cj; 
l  ^  6  ui  ^  if  ^  +  0  +  w<^>j<7r  +  u;. 

The  first  case  ir  +  uj<  fa  —  9  +  u  corresponds  to  the 
scenario  where  the  sensor  is  located  above  the  bullet  trajectory 
and  the  third  case  ^  +  9  +  ui  <  fa  <  tt  +  u>  corresponds 


’The  arctan  2  in  (2)  yields  —  -rr  <  <j>i  <  tt.  Thus  2tt  must  be  added  to  the 
arctan  2  result  to  obtain  a  positive  fa  if  fa  <  0. 


to  the  scenario  where  the  sensor  is  located  below  the  bullet 
trajectory  (as  shown  in  figure  1).  The  case  where  fa  =  n  +  u> 
corresponds  to  the  scenario  when  the  sensor  is  located  on  the 
bullet  trajectory  and  here  we  do  not  consider  such  a  scenario. 

If  fa  is  outside  the  bound  given  in  (4),  the  sensor  node  only 
receives  the  the  muzzle  blast  and  it  is  outside  the  FOV  of  the 
shockwave.  Under  the  assumptions  that  the  bullet  maintains 
a  constant  velocity  over  its  trajectory,  the  time  difference 
between  the  shockwave  and  the  muzzle  blast  can  be  written 
as  [2] 

n 

Tt  m  —  [1  -  cos\fa  -  tpi\\ ,  Vfa^<pi  (6) 

where  c  indicates  the  speed  of  sound.  Utilizing  (5),  the  bullet 
trajectory  angle,  lo,  can  be  obtained  from  the  shockwave  DOA 
angle.  Though  this  paper  assumes  that  the  bullet  speed  is 
constant  over  its  trajectory,  others  have  proposed  localization 
algorithms  [10],  [11]  that  employ  more  realistic  bullet  speed 
models  at  the  expense  of  computational  efficiency. 

III.  Data  Fusion  at  Sensor  Node  Level 
When  the  sensor  node  is  within  the  FOV  of  the  shockwave, 
the  three  available  measurements  are  the  two  DOA  angles  and 
the  time  difference  of  arrival  (TDOA)  between  the  muzzle 
blast  and  the  shockwave,  i.e., 

fa  =  h1(T:Siiij)  +  ri$  (7) 

fa  =h2(T,  Si,uj)  +  i)lp  (8) 

n  =  h3  (T,  Si,  fa  +  r?T  (9) 

where  /ii(-)  is  given  in  (2),  /i2 ( - )  is  given  in  (5),  and 
/z3 (-)  is  given  in  (6).  The  measurement  noise  is  assumed  to 
be  zero  mean  Gaussian  white  noise,  i.e.,  M{Q,a2fa, 

r/v  ~  A/'(0,ap  and  rjT  ~  Af(0,a2).  Let  T;  =  [fa  fa  fa] 
denotes  the  individual  sensor  level  estimates  on  the  target 
bearing,  range,  and  the  bullet  trajectory.  Data  fusion  at  the 
sensor  node  involves  calculating  these  individual  estimates 
based  on  the  three  sensor  measurements. 

Using  (5),  the  bullet  trajectory  angle,  ui,  can  be  obtained 
from  the  shockwave  DOA  measurements.  Thus,  the  observa¬ 
tions  on  the  trajectory  angle  can  be  written  as 

fa  =  u>i  +  r (10) 
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Now  the  likelihood  function,  p  (&>i\T,  Si,  u>),  can  be  written  as 
p(u)i\T,Si,u;)  =  Af{u,a2) 

From  (6),  the  range  can  be  written  in  terms  of  the  TDOA  as 

CT j 

n  =  Ti - rr - tt  (ID 

[1  -  cos  I  (pi  -  Pi\\ 

The  observation  of  r,  may  be  written  as 

„  CTi 

n  = 


1  —  cos 


■  (pi 


(12) 


Using  the  first-order  Taylor  series,  the  range  measurement  can 
be  approximated  as 

CTi 


[1  -  cos  I <j>i  -  Pi\] 


_ c _  CTi  sin\<j>j—(pi  | 

_  [  1  COS \(f>i  ^Pi  | ]  [l-COs|c/>i-^i|]: 

=  n  +  H(T,  Si,  U))rjr 


Vr 

T/rfxp 


where 


Vr  = 


Vr 

Thf>ip 


and 


&  Vtv  ~  V  (0,  <?l  +  <rl) 

zj(rr  Q.  ,  [ _ c _  CTi  sin|<ftj  —  (fi  | 

yl  j  Di,  UJ)  [1— cos\(j>i— 


~Vi\\  [1-C0S|^-^|7 

Now  the  likelihood  p  (fi|T,  Si,  to)  can  be  approximated  as 
p(h\T,  Si,u>)  «  Af  (n,<Jr(T,Si,u)j) 
where  the  variance  <J2r{T,  St,uj)  can  be  written  as 

7?  0 


a2(T,Si,co)  =  H(T,Si,co) 


0  al  +  al. 


(fi\ T,Si,J) 
p(fi\T,Si ,w)  «AT(MTi,ETj 


Thus,  the  likelihood  function  p 
mated  as 


HT(T,  Si,u>) 

(13) 

can  be  approxi- 

)  (14) 


where 

Mr*  = 


It  is  assumed  that  a  GPS  receiver  is  used  to  obtain  an  accurate 
positioning  on  each  sensor.  Thus,  the  position  observation  on 
the  sensors  are  given  as 


~al 

0 

0  ' 

Ti 

&  Sjt  — 

0 

*2AT,Si,u) 

0 

UJ 

_0 

0 

5,= 


'Six' 

Via 

+ 

Siy_ 

Viy_ 

(15) 


where  the  noise  terms  are  assumed  to  be  zero  mean  Gaussian 
white,  i.e.,  v ix  ~  Af(0,afx)  and  Viy  ~  J\f(0,afy).  Now  the 
GPS  measurement  likelihood  function  may  be  written  as 

t2 


p  (SilSj)  ~  AT 


0 


0 


=  Af  (vs, ,  ’Esi ) 


(16) 


Assumption  1.  Without  loss  of  generality,  it  can  be  assumed 
that  the  GPS  observations  on  sensor  position  are  independent 
of  target  location,  observations  on  target  location,  and  the 
projectile  trajectory  information,  i.e., 

p(Si|5i)  =p(Si\T,Si,Lo)  =p(Si\Ti,T,Si,uj) 

Base  on  assumption  1,  the  joint  probability 
p  (p'i,  Si\T,  Si,(jj^J  can  be  calculated  as 

P  ( Ti ,  Si\T ,  Si,u)  =  p  ( Si\Ti ,  T,  Si,u>)  p  (f \\T,  Si,u>) 

(17) 

Substituting  (14)  and  (16),  the  above  joint  likelihood  can  be 
written  as 

P  (Ti,  Si\T,  Si,u?j  «Ar(MSi,E5i)AA(MTi,ETi)  (18) 

Now  using  the  Bayes’  rule,  the  node  level  estimates  are 
given  as 

p(r,SiMTi,s^  = 

p(Ti,Si\T,Si,u)p(T,Si,u)  (19) 

fffp  (fi,  Si\T,  Si,ujp  (T,  Si,  u>)  dTdSiduj 

Note  that  the  denominator  in  (19)  indicates  the  normalization 
factor  and  since  no  a  priori  information  is  assumed  to  be 
known,  a  uniform  pdf  may  be  selected  for  p  ( T ,  Si,u>).  Since 
the  denominator  is  the  normalizing  term,  which  is  a  constant 
with  respect  to  T,  Si,  and  u>,  equation  (19)  can  be  written  as 

p(T,Si,u>\Ti,fy  «  ap  (fi,  Si\T,  Si,  ujj  (20) 

where  a  is  a  constant. 

Now  for  a  sensor  located  in  the  FOV  of  the  shockwave,  the 
target  location  can  be  estimated  as: 

TXi  =  Six  +  u  cos  (fi)  (21) 

Tyi  =  Siy  +  fi  sin(^j)  (22) 

When  the  sensor  is  located  outside  the  shockwave  FOV,  the 
only  estimate  would  be  the  bearing  angle.  After  individual 
estimates  are  obtained  at  the  sensor  node  level,  the  measured 
information  is  transmitted  to  a  central  node  where  it  is  fused 
to  obtain  a  more  accurate  estimate  of  shooter  location. 

IV.  Data  Fusion  at  the  Central  Node 

While  sensors  in  the  FOV  of  the  muzzle  blast  and  the  shock- 
wave  yield  a  range,  bearing,  and  trajectory  angle  estimates,  the 
gunfire  detection  systems  outside  the  FOV  of  the  shockwave 
yield  a  muzzle  blast  DOA.  Also,  GPS  measurements  are 
available  on  each  sensor  locations.  At  the  central  node,  this 
information  from  the  individual  sensor  nodes  is  fused  to  obtain 
an  accurate  estimate  of  the  shooter  location,  bullet  trajectory 
angle,  and  the  sensor  location. 
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Based  on  assumption  1,  the  joint  likelihood  function  associ¬ 
ated  with  each  sensor,  i.e.,  p  \  Ti ,  5*|T,  Si,uj^j,  can  be  written 
as 

P  (fi,  Si\T,  Sh  w)  =  p  (Silfi,  T,  Su  u)  p  (fi\T,  Su  w) 

Let  Si:n  =  {Sj,  S2,  —  ,  5n},  Ti:n  =  {Ti,  T2, . . . ,  Tn},  and 
Si:n  =  {Si,  S2,  ■  •  • ,  Sn},  where  n  indicates  the  number  of 
sensors.  Since  the  sensor  nodes  are  independent  of  each  other, 
the  joint  conditional  density  p  ^Ti:n,  Si:n|T,  can  be 

defined  as 

n 

P  (T1:n,  Si:n|T,  S1:n,  ojJ  =Hp(fz,Si\T,Si,^  (23) 


Here  we  consider  the  maximum  likelihood  approach  to  obtain 
the  fused  estimates.  In  the  maximum  likelihood  estimation  ap¬ 
proach  considered  here,  an  estimate  of  the  sensor  locations,  the 
shooter  location,  and  the  bullet  trajectory  angle  are  obtained 
so  that  the  joint  log-likelihood  function  is  maximized,  i.e.. 


max 

T,Sl:„,0! 


In  jp  (T1:n,S1:n|T,  Si:„,cu)  j 

n 

=>■  ^max  ^2  In  {p  (fi,  S, i\T ,  si,  } 


Based  on  the  results  given  in  the  previous  section,  the  criteria 
for  the  maximum  likelihood  estimation  can  be  written  as 


n 

max  Iln  W  (VTi,  StJ)  +  In  {A f  (nSi,  S5J}] 

T,Sl:n,CO  ^ 

(25) 


Note  that  the  densities  A /(/xt*,  EtJ  and  A/”  {fJ-Si,  EsJ  may 
be  written  as 


Af(,‘T"ET-) = 7raexp{ ~  5 

(fi  -  l-l 7\ )  J 


(26) 


where 


<t>i 

arctan  2  {Tv  -  Siy ,  Tx  -  S ^ ) 

IJ-Ti  = 

n 

= 

V(; TX-SiJ2  +  (Ty-Siy )2 

U) 

CO 

& 


S71.  = 


0 


0 


0  a2r{T,Si,uj)  0 


0 


0 


and 


Af("s'Ss-) = 7raexp{  ■  5  (*•-«'/ 

^s!  (&i  -  MSi)  | 


(27) 


where 


'Si  1 

[of  0  1 

MSi  = 

kJ 

&  E  Ti  — 

'lx 

0  a2 

L  Ly  J 

After  substituting  (26)  and  (27)  into  (25),  the  maximum 
likelihood  criteria  may  be  written  as 


n 

min  > 

T,Si-n,co  ^ 

1 

1  1 

to  | 

(Ti  —  flTi 

)Tk 

(f 

~  T  t,  ) 

1 

/ft 

)Tk 

\ 

2 

ySi  -  usi 

(sz 

-  VSi) 

In  {V|27rETi|}+ In  {^2^1} 

(28) 


Note  that  the  term,  I11 j  \/j27rEj’i  ||,  in  above  equation  is 
present  due  to  the  fact  that  E-/;.  is  a  function  of  T,  S,  and 
uj.  The  last  term,  In  |  \/|27rEsj|  can  be  ignored  since  E 
is  a  known  constant  matrix.  Since  E^  is  assumed  to  be  a 
diagonal  matrix,  (28)  can  be  rewritten  as 


min 

T,  Sl:„,W 


E 


1  l 


ln(eoy)  +  ^  (fi 


where  e  is  defined  as 


e  =  (2n)3/2a(j,av, 

Apart  from  the  initial  term,  lnjeoy),  the  optimization  problem 
given  in  (29)  is  similar  to  that  used  in  the  weighted  nonlinear 
least-squares.  Thus,  the  maximum  likelihood  approach  pre¬ 
sented  here  is  similar  to  the  weighted  nonlinear  least-squares 
estimation. 


V.  Nonlinear  Least  Squares 

There  exist  no  closed  form  solution  to  the  nonlinear  least- 
squares  optimization  problem  given  in  (29)  and  therefore  a  nu¬ 
merical  approach  needs  to  be  used.  A  few  common  approaches 
to  solving  the  nonlinear  least-squares  problem  include  the 
Gauss-Newton  method,  Nelder-Mead  simplex  method,  and 
Marquardt  method  [12].  Almost  all  these  approaches  are 
iterative  methods  that  require  an  initial  approximation  to  the 
unknown  parameters  and  provide  successively  better  approxi¬ 
mations.  The  iterative  process  is  repeated  until  the  parameters 
do  not  change  to  within  specified  limits. 

This  section  provides  the  Gauss-Newton  method  for  solving 
the  nonlinear  least  squares  problem  given  in  (29).  The  main 
advantage  of  the  Gauss-Newton  method  is  that  it  exhibits 
a  “quadratic  convergence,”  which,  simply  put,  means  that 
the  uncertainty  in  the  parameters  after  p  +  1  iterations  is 
proportional  to  the  square  of  the  uncertainty  after  p  iterations. 
Once  these  uncertainties  begin  to  get  small,  they  decrease  quite 
rapidly.  An  additional  advantage  of  the  Gauss-Newton  method 
is  that  it  only  requires  calculating  the  first-order  derivatives. 
The  major  problem  with  the  Gauss-Newton  method  is  that  it 
sometimes  diverges  if  the  initial  approximation  is  too  far  from 
truth. 
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In  order  to  simplify  the  formulation,  we  treat  Ey  as  a 
known  constant  matrix.  Thus,  (29)  can  be  rewritten  as 


1  rr 

min  J  =  -  Ay  T IV  Ay 


(30) 


■fl 

~  Mr, 
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0 

0 

0 

0  ' 

-  MSi 
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0 

0 

0 

,W  = 

0 

0 

0 

0 

Tn 

- 

0 

0 

0 

0 

s„ 

-  /Ls„_ 

.  0 

0 

0 

o" 

y-t 

where 


Ay  = 


Let  x  =  [Tx  Ty  ui  Sla  Sly  ...  Snx  S„y]  de¬ 
note  the  parameters  to  be  estimated  and  let  y  = 
[f^ (x)  ...  fzT  (x)  ...  f,J (x)]  denote  the  measure¬ 

ments.  Here  f;(x)  =  \tj>i  ri  <o  Six  5iB]T  is  defined  as 

arctan  2  (Ty  -  Siy ,  Tx  -  Six) 

y/(Tx  -  si)2  +  (Ty  -sly)2 

UJ 

Si . 

Also  let  y  =  [yf , . . .  ,yf , . . .  ,y£]T  where  y* 

.  .if 

T'f  Sj  ■  Now  Ay  is  defined  as 
Ay  =  y  -  y 

Let  the  current  estimates  estimates  of  x  be  denoted  as 


f i  (x)  = 


(31) 


xc  =  [Tcx  T- 

Define 


wc  siv 

Ax  =  x  -  xc 


If  the  components  of  Ax  are  sufficiently  small,  then  using  the 
first-order  Taylor  series  approximation,  we  have 


Figure  3.  Gauss-Newton  Algorithm. 


Now  using  the  current  estimates,  F,  W,  and  Ayc  are  calcu¬ 
lated.  Then,  Ax  estimate  for  the  next  iteration  is  calculated 
from  (35)  and  this  process  is  repeated  until  Ax  converges  to 
a  prescribed  small  value.  A  schematic  representation  of  the 
Gauss-Newton  algorithm  is  presented  in  figure  3. 

VI.  Results 


f  (x)  «  f  (xc)  -I-  F Ax 


(32) 


where 


'=i 


Now  the  measurement  residual  Ay  can  be  linearly  approxi¬ 
mated  as 


This  section  presents  numerical  simulations  to  assess  the  lo¬ 
calization  improvement  due  to  the  proposed  fusion  algorithm. 
Here  we  consider  two  separate  simulation  scenarios,  for  both 
scenarios,  we  assume  that  there  are  five  sensor  nodes  located 
at 


127  20  90  136  182 
107  22  0  68  59 


Ay  «  y  -  f  (xc)  -  F Ax  =  Ayc  -  F Ax  (33) 

where  Ayc  =  y  —  f  (xc).  Substituting  (33)  in  (30)  yields 

J  «  i  ^ Ayc  -  FAy}j  W  (^Ayc  -  LAxj  (34) 

The  Ax  that  minimizes  the  above  cost  function  can  be  written 
as 

Ax=  {FTWF)~1  FTWAyc  (35) 

After  obtaining  Ax,  the  current  estimates  are  redefined  as 

xc  =  Ax  +  xc  (36) 


For  simulation  purposes,  we  assume  a  constant  velocity  model 
for  the  bullet.  Thus,  the  Mach  number  is  selected  to  be  m  =  2 
and  the  speed  of  sound  is  selected  to  be  c  =  342  m/ sec.  For 
both  scenarios,  the  measurement  noise  models  are  selected  as 
=  <Ay  =  5 m,  <70  =  av  =  4°,  and  aT  =  1  msec.  Since 
there  exist  several  approaches  to  solve  the  nonlinear  least- 
squares  problem,  two  different  methods  are  used  to  obtain 
solutions  for  both  simulation  scenarios.  In  the  first  method, 
the  optimization  problem  is  solved  using  the  Gauss-Newton 
method  [12]  presented  in  the  previous  section.  The  second 
approach  uses  the  Nelder-Simplex  algorithm  [13],  i.e.,  the 
fminsearch  function  in  Matlab. 
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A.  Simulation  I 

For  the  first  simulation,  the  shooter  is  assumed  to  be  located 
at  T  =  [50  50]  and  we  select  the  bullet  trajectory  to  be 

u>  =  30°.  Figure  4  shows  the  first  simulation  scenario.  Due  to 
the  sensor  locations,  the  second  and  the  third  sensors  do  not 
receive  the  shockwave. 


a 


-  -  -  Bullet  Trajectory 

Shooter  Location 

•  Sensor  Locations 

-  Shock  FOV 

-50 - 1 - 1 - 1 - 1 - 1 - 

0  50  100  150  200  250  300 

ai-axis 

Figure  4.  Simulation  I  Scenario 


In  order  to  evaluate  the  system  performance,  a  Monte  Carlo 
simulation  is  conducted  for  both  the  Gauss-Newton  method 
and  the  simplex  algorithm.  The  mean  shooter  locations  and 
the  associated  error  ellipses  obtained  from  the  Monte  Carlo 
simulations  using  the  Gauss-Newton  method  are  given  in 
figure  5.  A  separate  plot  is  not  provided  for  the  results  obtained 
using  the  simplex  algorithm  since  they  are  very  similar  to  that 
obtained  for  the  Gauss-Newton  method.  Figure  5  indicates  that 
the  sensor  five  performs  the  worst  out  of  the  three  sensors  in 
the  shockwave  FOV.  Figure  5  also  indicates  that  the  fused 
estimate  is  superior  to  the  individual  sensor  estimates,  and 
the  uncertainty  associated  with  the  fused  estimates  is  much 
less  than  the  uncertainty  associated  with  the  individual  sensor 
estimates. 


60 
55  - 

■S3  so  - 

»«- 
40  - 
35  - 
30  - 


25  — 
-40 


★  Truth 

#  Fused 

■ - Fused  Err 

*  Sensor  1  Est 

-  -  ■  Sensor  1  Err 
$  Sensor  4  Est 

—  Sensor  4  Err 
$  Sensor  5  Est 

. Sensor  5  Err 


40  . 
x-axis 


individual  sensors  and  the  fusion  algorithms  over  the  Monte 
Carlo  run.  The  “average”  estimate  presented  in  table  I  indicates 
the  estimate  obtained  by  simply  averaging  the  individual  target 
estimate  from  sensors  one,  four,  and  five.  Table  I  also  contains 
the  root  mean  square  error  (RMSE)  associated  with  each 
estimate.  Based  on  the  RMSE  presented  in  table  I,  one  could 
conclude  that  that  fused  estimates  outperform  the  individual 
sensors  and  the  simple  average  estimate. 

Table  I 

Simulation!:  Shooter  Location 


Tx  (m) 

Tv  (m) 

RMSE  (m) 

Truth 

50 

50 

- 

Sensor  1 

48.3513 

47.2948 

23.2870 

Sensor  2 

- 

- 

- 

Sensor  3 

- 

- 

- 

Sensor  4 

42.9248 

50.2141 

31.1132 

Sensor  5 

37.1197 

52.0782 

65.6542 

Average 

42.7986 

49.8623 

25.9660 

Gauss-Newton 

49.9066 

49.9134 

6.8639 

Nedler-Simplex 

50.0493 

50.0588 

6.9972 

Table  II  contains  the  mean  bullet  trajectory  angle  estimate 
obtained  from  the  individual  sensors  and  the  fusion  algorithms 
over  the  Monte  Carlo  run.  Table  II  also  contains  the  RMSE 
associated  with  each  trajectory  angle  estimate.  Note  that  the 
fused  trajectory  estimate  is  simply  the  average  of  the  individual 
sensor  estimates  due  to  the  way  in  which  uj  appears  in  (31). 


Figure  6.  Simulation  I:  Localization  error  histogram  for  average  estimates 

Figure  6  presents  the  shooter  localization  error  histogram  for 
the  average  estimate,  i.e.,  the  estimate  obtained  by  simply  av¬ 
eraging  the  individual  target  estimates  from  sensors  one,  four, 
and  five.  Figure  7  presents  the  localization  error  histogram  for 
the  fused  estimate  obtained  for  the  Gauss-Newton  method. 

Table  III  contains  RMSE  associated  with  the  sensor  location 
estimates.  Interestingly,  the  fusion  algorithm  was  able  to 
improve  the  sensor  location  accuracy  by  reducing  the  GPS 
uncertainties.  Based  on  the  RMSE  presented  in  tables  I,II,  and 
III,  one  could  conclude  that  that  fused  estimates  outperform 
the  individual  sensors. 


Figure  5.  Simulation  I:  Mean  Results  from  Monte  Carlo  Runs 
Table  I  contains  the  mean  shooter  location  estimate  of  the 


B.  Simulation  II 

For  the  second  simulation,  the  shooter  is  assumed  to  be 
located  at  T  =  [150  —  50]T  and  we  select  the  bullet 
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Figure  7.  Simulation  I:  Localization  error  histogram  for  fused  estimates 
Table  II 

Simulation!:  Bullet  Trajectory 


u  (deg) 

RMSE  (deg) 

Truth 

30 

- 

Sensor  1 

30.0641 

3.9690 

Sensor  2 

- 

- 

Sensor  3 

- 

- 

Sensor  4 

30.3402 

3.9970 

Sensor  5 

29.9591 

3.9029 

Average 

30.1211 

2.2128 

Gauss-Newton 

30.1211 

2.2128 

Nedler-Simplex 

30.1999 

2.4674 

Table  III 

Simulation  I:  Sensor  Location  RMSE 


GPS  (m) 

Gauss-Newton  (m) 

Nedler-Simplex  (m) 

Sensor  1 

7.0215 

6.5453 

6.5938 

Sensor  2 

7.0002 

6.3195 

6.3530 

Sensor  3 

7.0028 

6.6513 

6.6770 

Sensor  4 

7.1509 

6.5259 

6.6201 

Sensor  5 

7.0223 

6.7883 

6.8731 

both  sensor  two  and  sensor  three  are  of  similar  accuracy  since 
they  are  equal  distance  from  the  bullet  trajectory.  Figure  9 
also  indicates  that  the  fused  estimate  is  superior  to  individual 
sensor  estimates  and  the  uncertainty  associated  with  the  fused 
estimates  is  much  less  than  the  uncertainty  associated  with  the 
individual  sensor  estimates. 


Figure  9.  Simulation  II:  Mean  Results  from  Monte  Carlo  Runs 

Table  IV  contains  the  mean  shooter  location  estimate  of  the 
individual  sensors  and  the  fusion  algorithms  over  the  Monte 
Carlo  runs.  Table  IV  also  contains  the  RMSE  associated  with 
each  estimate.  Based  on  the  RMSE  presented  in  table  IV, 
one  could  conclude  that  that  fused  estimates  outperform  the 
individual  sensors  and  the  simple  average  estimate. 


trajectory  to  be  u>  =  170°.  Figure  8  shows  scenario  for  the 
second  simulation.  As  shown  in  figure  8,  only  the  second  and 
the  third  sensors  are  in  the  FOV  of  the  shockwave. 
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Figure  8.  Simulation  II  Scenario 


The  mean  shooter  locations  and  the  associated  error  ellipses 
obtained  from  the  Monte  Carlo  simulation  using  the  Gauss- 
Newton  method  are  given  in  figure  9.  Figure  9  indicates  that 


Table  IV 

Simulation  II:  Shooter  Location 


Tx  (m) 

Tv  (m) 

RMSE  (m) 

Truth 

150 

-50 

- 

Sensor  1 

- 

- 

- 

Sensor  2 

156.5714 

-56.8403 

49.9558 

Sensor  3 

156.3943 

-57.4820 

43.3278 

Sensor  4 

- 

- 

- 

Sensor  5 

- 

- 

- 

Average 

156.4828 

-57.1611 

34.2454 

Gauss-Newton 

150.1808 

-50.6793 

10.3593 

Nedler-Simplex 

150.0338 

-50.5054 

10.7630 

Figure  10  presents  the  shooter  localization  error  histogram 
for  the  average  estimate  and  figure  1 1  presents  the  localization 
error  histogram  for  the  fused  estimate  obtained  for  the  Gauss- 
Newton  method. 

Table  V  contains  the  mean  bullet  trajectory  angle  estimate 
and  associated  RMSE  obtained  from  the  individual  sensors 
and  the  fusion  algorithms  over  the  Monte  Carlo  runs.  Finally, 
table  VI  contains  RMSE  associated  with  the  sensor  location 
estimates  for  simulation  two.  Note  that  the  fusion  algorithm 
was  able  to  improve  the  sensor  location  accuracy  by  reducing 
the  GPS  uncertainties.  The  RMSE  presented  in  tables  IV, V, 
and  VI,  indicate  that  fused  estimates  outperform  the  individual 
sensors  for  the  second  simulation. 
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Figure  10.  Simulation  II:  Localization  error  histogram  for  average  estimates 


Figure  1 1 .  Simulation  II:  Localization  error  histogram  for  fused  estimates 


Table  V 

Simulation II:  Bullet  Trajectory 


U)  (deg) 

RMSE  (deg) 

Truth 

170 

- 

Sensor  1 

- 

- 

Sensor  2 

169.8829 

4.1206 

Sensor  3 

169.9445 

4.0493 

Sensor  4 

- 

- 

Sensor  5 

- 

- 

Average 

169.9137 

2.8139 

Gauss-Newton 

169.9137 

2.8139 

Nedler-Simplex 

169.9098 

3.0287 

Table  VI 

Simulation  II:  Sensor  Location  RMSE 


GPS  (m) 

Gauss-Newton  (m) 

Nedler-Simplex  (m) 

Sensor  1 

6.9913 

6.7276 

6.7284 

Sensor  2 

6.9897 

6.7586 

6.7848 

Sensor  3 

6.9484 

6.6360 

6.6923 

Sensor  4 

7.0946 

6.7296 

6.7282 

Sensor  5 

7.0078 

6.7280 

6.7317 

VII.  Final  Remarks 

The  shooter  localization  problem  using  a  network  of  soldier- 
worn  gunfire  detection  systems  is  considered  here.  This  paper 
presents  a  fusion  algorithm  that  utilizes  the  benefits  of  the 
sensor  network  layout  of  all  the  sensors  within  a  small 
combat  unit  to  help  refine  the  shooter  localization  accuracy. 
The  individual  gunfire  detection  systems  considered  here  are 
composed  of  a  passive  array  of  microphones  that  is  able  to 


localize  a  gunfire  event  by  measuring  the  direction  of  arrival 
for  both  the  muzzle  blast  and  the  shockwave.  After  detecting 
a  gunfire,  the  individual  sensors  report  their  solution  along 
with  their  GPS  positions  to  a  central  fusion  node.  At  the 
central  node,  the  individual  solutions  are  fused  along  with 
the  GPS  locations  on  the  soldiers  to  yield  a  highly  accurate 
geo-rectified  solution.  Numerical  results  given  here  indicate 
that  the  fused  estimates  are  more  accurate  than  the  individual 
localization  results.  Future  work  include  further  analyzing  the 
linearization  issues  associated  with  the  maximum  likelihood 
approach  and  developing  a  mathematically  rigorous  method 
to  quantify  the  uncertainties  associated  with  the  maximum 
likelihood  estimates. 
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